In recent years the R spatial community has had a big impact on big data, especially with the processing of spatial data, both in vector format (points, lines, polygons, etc.) and raster format (satellite images, drones, etc.).

Nowadays, there are many packages within R spatially to work on a key aspect within spatial analysis, from the evaluation of a Moran test to cluster identification, geographic weighted regression, pre and post processing of satellite images or drone images, among others; there is also a great potential within spatial data visualizations, from simple static visualizations, to dynamic and interactive ones, including 3D visualizations.

However, there is a huge gap for processing large volumes of data without high computational costs, and it is precisely in this aspect that the rgee package aims to work and be part of this ecosystem to process large spatial datasets in an accessible and user-friendly way hand in hand with the integration of the spatial universe of R packages.

Mapping the past temperature and the current vegetation index

For this example, we use the ERA5-Land dataset, a reanalysis dataset providing a consistent view of the evolution of land variables over several decades at an enhanced resolution compared to ERA5, more information click here. together MOD13A1 V6 product provides a Vegetation Index, more information click here.

Information:

- For this demo you need rgee, tmap, raster, viridis,cptcity and mapview packages previously installed.

3.1 Requeriments

library(rgee)
library(tmap)
library(raster)
library(viridis)
library(mapview)
library(cptcity)
ee_Initialize() # Initialize gee inside R with registered account

── rgee 1.1.5 ────────────────────────────────────────────────────── earthengine-api 0.1.334 ── 
 ✔ user: not_defined
 ✔ Initializing Google Earth Engine:  DONE!
 ✔ Earth Engine account: users/ambarja 
───────────────────────────────────────────────────────────────────────────────────────────────

3.2 Processing data with rgee

r_temp <- ee$ImageCollection$Dataset$ECMWF_ERA5_LAND_MONTHLY
ee_temp <- r_temp$
  select("temperature_2m")$
  filterDate("1992-01-01","1992-12-31")$
  mean()$
  subtract(273.15)

3.3 Simple visualization

viz_temp <- list(
  min = 5,
  max = 25,
  palette = magma(n = 100)
  )

temp <- Map$addLayer(ee_temp, visParams = viz_temp) + 
  Map$addLegend(visParams = viz_temp,name = "Tem[C°]",position = "bottomleft")

temp

3.4 Define box to study area

geometry <- ee$Geometry$Rectangle(
  coords = c(-180,-90,180,90),
  proj = "EPSG:4326",
  geodesic = FALSE
)

3.5 Identifying min and max values of temperature

(minmax <- ee_temp$reduceRegion(
  reducer = ee$Reducer$minMax(),
  geometry = geometry,
  scale = 500*1000)$
  getInfo())

$temperature_2m_max
[1] 29.68718
$temperature_2m_min
[1] -53.14849

3.6 Generate a thumbnail image in a rasterlayer or stars format

world_temp <- ee_as_thumbnail(
  image = ee_temp,     # Image to export
  region = geometry,   # Region of interest
  dimensions = 1024,   # output dimension
  raster = TRUE,       # If FALSE returns a stars object. FALSE by default
  vizparams = list(min = -53.14849, max = 29.68718)
)

3.7 Processing data for mapping with tmap

# Define a projection (EPSG:4326 - WGS 84)
crs(world_temp) <- 4326
world_temp[] <- scales::rescale(
  getValues(world_temp),
  c(minmax$temperature_2m_min, minmax$temperature_2m_max)
) 
# Load world vector (available after load tmap)
data("World") 
# Define a Robin type projection 
robin <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
# Remove the temperature value of sea
map <-  world_temp %>% 
  mask(World)

3.8 Exploration raster data

tm_shape(map) +
  tm_raster(palette = viridis(n = 100),style = "cont")

3.9 Elaboration of the Temperature Map

temp_map <- map %>% 
  tm_shape(projection = robin,raster.warp = FALSE) +
  tm_raster(
    title = "Temperature (°C)",
    palette = magma(n = 100),
    style = "cont"
  ) + 
  tm_shape(shp = World) + 
  tm_borders(col = "black",lwd = 0.7) + 
  tm_graticules(alpha = 0.8, lwd = 0.5, labels.size = 0.5)
temp_map

3.10 Export map in png format

tmap_save(temp_map,"temp_map.png",width = 6,height = 3)

3.11 Mapping of the current vegetation index

# Processing data with rgee
r_ndvi <- ee$ImageCollection$Dataset$MODIS_006_MYD13Q1
ee_ndvi <- r_ndvi$
  select("NDVI")$
  filterDate("2022-01-01","2022-01-31")$
  max()$
  multiply(0.0001)

3.12 Simple visualization

viz_ndvi <- list(
  min = -0.5 ,
  max = 0.87,
  palette = cpt(pal = "grass_ndvi")
  )
ndvi <- Map$addLayer(ee_ndvi,visParams = viz_ndvi) + 
  Map$addLegend(visParams = viz_ndvi,name = "NDVI")
ndvi